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ABSTRACT 

We study the inhomogeneous reionization in a critical density CDM universe due to stellar sources, 
including Population III objects. The spatial distribution of the sources is obtained from high resolution 
numerical N-body simulations. We calculate the source properties taking into account a self-consistent 
treatment of both radiative {i.e. ionizing and H2 -photodissociating photons) and stellar {i.e. SN explo- 
sions) feedbacks regulated by massive stars. This allows us to describe the topology of the ionized and 
dissociated regions at various cosmic epochs and derive the evolution of H, He, and H2 filling factors, 
soft UV background, cosmic star formation rate and the final fate of ionizing objects. The main results 
are: (i) galaxies reionize the IGM by zfulO (with some uncertainty related to the gas clumping factor), 
whereas H2 is completely dissociated already by z«25; (ii) reionization is mostly due to the relatively 
massive objects which collapse via H line cooling, while objects whose formation relies on H2 cooling 
alone are insufficient to this aim; (iii) the diffuse soft UV background is the major source of radiative 
feedback effects for z < 15; at higher z direct flux from neighboring objects dominates; (iv) the match 
of the calculated cosmic star formation history with the one observed at lower redshifts suggests that 
the conversion efficiency of baryons into stars is «!%; (v) we find that a very large population of dark 
objects which failed to form stars is present by z«8. We discuss and compare our results with similar 
previous studies. 

Subject headings: galaxies: formation - cosmology: theory 



1. INTRODUCTION 

At zwllOO the intergalactic medium (IGM) is expected 
to recombine and remain neutral until the first sources of 
ionizing radiation form and reionize it. The application of 
the Gunn-Peterson (1965) test to QSOs absorption spec- 
tra suggests that the HI reionization is complete by z«5. 
Recently, Lya emitters have been seen up to a redshift 
of z=6.68 (Chen, Lanzetta & Pascarelle 1999). This in- 
dicates that the IGM was already ionized by then, but 
the epoch of complete reionization has yet to be firmly 
established. Several authors (Shapiro, Giroux & Babul 
1994; Madau, Haardt & Rees 1998 and references therein) 
have claimed that the known population of quasars (show- 
ing a pronounced cutoff a,t z ^ 3 [Warren, Hewett & Os- 
mer 1994; Schmidt, Schneider & Gunn 1995; Shaver et 
al. 1996]) and galaxies provides «10 times fewer ionizing 
photons than are necessary to keep the observed IGM ion- 
ization level. Thus, additional sources of ionizing photons 
are required at high redshift. Although models have been 
proposed in which the UV photons responsible for the IGM 
reionization may be emitted by a cosmological distribution 
of decaying dark matter particles, such as neutrinos (see 
for example Scott, Rees & Sciama 1991), the most promis- 
ing sources are early galaxies and quasars. 

Recent observational evidence suggest the existence of 
an early population of pregalactic objects which could have 
contributed to the reionization and metal enrichment of 
the IGM. Metals have been clearly detected in Lya forest 
clouds with column densities low enough to be identified 



with a truly diffuse IGM (Cowie et al. 1995; Tytler et 
al. 1995; Lu et al. 1998; Cowie k Songaila 1998), although 
discrepant metallicities have been inferred. Naively, one 
can estimate that the amount of heavy elements associ- 
ated with the number of photons required to ionize ev- 
ery baryon in the universe a few times corresponds to an 
IGM metallicity Z « 3 x 10'^ Zq. This figure is roughly 
consistent with the value obtained by Cowie & Songaila 
(1998), but ten times larger than the one derived by Lu 
et al. (1998). The reasons for this discrepancy are not yet 
clear. Whatever the correct level of metallicity is, these 
heavy elements must have been produced by objects host- 
ing the sites of the earliest star formation activity in the 
universe. The presence of a soft, stellar component of 
the UV background deduced from studies of the [Si/C] 
abundance ratios in low column density absorption sys- 
tems (Savaglio et al. 1997; Giroux & ShuU 1997), also 
supports this view. The question remains concerning the 
transport mechanism from the production regions, which 
are presumably associated with high peaks of the density 
fluctuation field, into the very diffuse medium probed by 
absorption line experiments. It is currently unclear if the 
mechanical energy input associated with the same massive 
stars producing the metals {i.e. supernovae) is sufficient to 
drive the metals far enough from their production site or 
if additional mechanisms need be invoked. Although re- 
cently Gnedin (1998), using high resolution cosmological 
simulations which include a coarse modelling of a two- 
phase interstellar medium, has suggested that merging is 
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the dominant mechanism for transporting heavy elements 
from primeval galaxies into the IGM, direct ejection of 
the metal enriched interstellar gas by supernovae or stel- 
lar winds may still play an important role. 

A contribution to the ionizing flux at high rcdshift from 
an early population of quasars may provide the hard com- 
ponent required for Hell reionization. High resolution 
spectra of a quasar at zf»3 have shown large fluctuations 
of the HI to HeH optical depth ratio (Reimers et al. 1997), 
interpreted as evidence for a patchy HeH ionization in the 
IGM. However, it is unclear if these fluctuations could be 
rather caused by statistical fluctuations of the IGM den- 
sity or of the ionizing background flux (Miralda-Escude 
1998; Miralda-Escude, Haehnelt fc Rees 1998). Analytical 
results from Haiman & Loeb (1998b) show that Hell reion- 
ization due to sources with quasar-like spectra should oc- 
cur at approximately the same redshift as hydrogen reion- 
ization; if the predicted Hell reionization redshift will be 
confirmed by future measurements, this would severely 
constrain the number of high redshift quasars predicted 
by theory. 

The study of the IGM reionization due to these primeval 

sources has been tackled by several authors, both via 
semi-analytical (Shapiro 1986; Fukugita & Kawasaki 1994; 
Tegmark, Silk & Blanchard 1994; Haiman & Loeb 1998b; 
Miralda-Escude, Haehnelt & Rees 1998; Valageas & Silk 
1999, VS) and numerical (Gnedin & Ostriker 1997, GO; 
Norman, Paschos & Abel 1998; Razoumov & Scott 1998) 
approaches. One of the main difficulties of the problem 
is the treatment of cosmological radiative transfer. Only 
recently have numerical approaches have been attempted 
(Abel, Norman & Madau 1998; Norman, Paschos & Abel 
1998), and these still require severe approximations to the 
radiative transfer equation. Clearly, a numerical approach 
off'ers the advantage of a more detailed treatment of fea- 
tures such as the IGM dumpiness (GO; Shapiro, Raga & 
Mellema 1998), or the source spatial distribution. Con- 
versely, semi-analytical approaches are more; flexible and 
allow a more thorough exploration of the parameter space. 

A proper treatment of the reionization problem must 
take into account the interplay of galaxy formation with 
the reionization process itself. Basically, two types of feed- 
back, radiative and stellar, can be at work. In the stan- 
dard cosmological hierarchical scenario for structure for- 
mation, the objects which form first are predicted to have 
masses corresponding to virial temperatures Tyir < 10** 
K. Once the gas has virialized in the potential wells of 
dark matter halos, additional cooling is required to further 
collapse the gas and form stars. For a gas of primordial 
composition at such low temperatures the main coolant is 
molecular hydrogen (Peebles & Dickc 1968; Shapiro 1992; 
Haiman, Rees & Loeb 1996; Abel et al. 1997a; Tegmark 
et al. 1997; Ferrara 1998). We define Pop III objects as 
those for which H2 cooling is required for collapse. Af- 
ter a H2 molecule gets rotationally or vibrationally ex- 
cited through a collision with an H atom or another H2 
molecule, a radiative de-excitation leads to cooling of the 
gas. This mechanism has been studied in great detail by 
several authors (Lepp & ShuU 1984; Hollenbach & McKee 
1989; Martin, Schwarz & Mandy 1996; Galh & Palla 1998) 
and has been applied to the study of the formation of pri- 
mordial small mass objects (Haiman, Thoul & Loeb 1996; 
Tegmark et al. 1997; Abel et al. 1997b). Primordial H2 



forms with a fractional abundance of w 10^^ at redshifts 
400 via the H^ formation channel. At redshifts ^ 110, 
when the Cosmic Microwave Background radiation (CMB) 
intensity becomes weak enough to allow for significant for- 
mation of H~ ions, a primordial fraction of /^f^ 2 x 10~^ 
(Shapiro 1992; Anninos & Norman 1996) is produced for 
model universes that satisfy the standard primordial nu- 
cleosynthesis constrain Q^h^ = 0.0125 (Copi, Schramm & 
Turner 1995), where fib is the baryon density parameter 
and Hq ~ lOOh km s~^ Mpc^^ is the Hubble constant. 
This primordial fraction is usually lower than the one re- 
quired for the formation of Pop HI objects, but during the 
collapse phase, the molecular hydrogen content can reach 
high enough values to trigger star formation. On the other 
hand, objects with virial temperatures (or masses) above 
that required for the hydrogen Lya line cooling to be ef- 
ficient, do not rely on H2 cooling to ignite internal star 
formation. As the first stars form, their photons in the en- 
ergy range 11.26-13.6 eV are able to penetrate the gas and 
photodissociate H2 molecules both in the IGM and in the 
nearest collapsing structures, if they can propagate that 
far from their source. Thus, the existence of a UV flux 
below the Lyman limit due to primordial objects, capable 
of dissociating the H2, could strongly influence subsequent 
small structure formation. Haiman, Rees & Loeb (1997, 
HRL), for example, have argued that Pop III objects could 
depress the H2 abundance in neighbor collapsing clouds, 
due to their UV photodissociating radiation, thus inhibit- 
ing subsequent formation of small mass structures. On 
the other hand, Ciardi, Ferrara & Abel (1999, CFA) have 
shown that the "soft-UV background" (SUVB) produced 
by Pop His is well below the threshold required for nega- 
tive feedback to be effective earlier than zss20. In princi- 
ple, the collapse of larger mass objects can also be influ- 
enced by an ionizing background, as gas in halos with a 
circular velocity lower than the sound speed of ionized gas 
may be prevented from collapsing due to pressure support 
in the gravitational potential. Thoul & Weinberg (1996) 
(see also Babul & Rees 1992) have however shown that the 
collapse is only delayed by this process. We will refer to 
this complex network of processes as "radiative feedback" . 

The other feedback mechanism is related to massive 
stars and for this reason we will call it "stellar feedback" . 
Once star formation begins in the central regions of col- 
lapsed objects it may strongly influence their evolution 
via the effects of mass and energy deposition due to mas- 
sive stars through winds and supernova explosions. These 
processes may induce two essentially different phenomena. 
Low mass objects are characterized by shallow potential 
wells in which the baryons are only loosely bound and a 
relatively small energy injection may be sufficient to expel 
the entire gas content back into the IGM, i.e. a blowaway, 
thus quenching star formation. Larger objects, may in- 
stead be able to at least partially retain their baryons, 
although a substantial fraction of the latter are lost in an 
outflow, i.e. a blowout (Ciardi & Ferrara 1997; Mac Low 
& Ferrara 1999, MF; Ferrara & Tolstoy 1999, FT). How- 
ever, even in this case the outflow induces a decrease of 
the star formation rate due to the global heating and loss 
of the galactic ISM. Omukai & Nishi (1999) have pointed 
out that the ionizing radiation of the flrst stars formed in 
Pop III objects can also produce an abrupt interruption of 
the star formation by dissociating the internal H2 content. 
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The effect of this feedback is very similar to tlie blowaway. 
as both processes are regulated by massive stars. 

Our aim here is to study and describe in detail the reion- 
ization process of the IGM. As the ionizing sources are not 
spatially homogeneously distributed, the pattern of the 
ionized regions is not uniform, i.e. inhomogeneous reion- 
ization takes place. Thus, a crucial ingredient of such cal- 
culations is the realistic description of the ionized region 
topology in the universe. This can be achieved only by nu- 
merical simulations which track the formation and merging 
of dark matter halos. For our purposes, these simulations 
must reach the highest possible mass resolution in order 
to unambiguously identify the very small Pop III objects 
initiating the reionization process. The second fundamen- 
tal ingredient is a proper treatment of radiative and stel- 
lar feedback processes on which a great deal of previous 
work, either by ourselves or other groups, has been accu- 
mulating ranging from semi- analytical studies, to numer- 
ical models and hydrodynamical simulations, particularly 
with regard to stellar feedback. Improving on such results 
we have been able to effectively model and include in the 
computations a wide number of effects governed by ioniz- 
ing/dissociating radiation and massive star energy injec- 
tion which will be presented and discussed in detail. The 
delicate link between the properties of dark matter halos, 
which can be reliably derived from N-body simulations, 
and those describing their baryonic counterparts is repre- 
sented by the prescription used to model star formation. 
The advantages that the proposed approach brings along 
should become evident by reading the paper: it allows to 
derive self-consistent inhomogeneous reionization models 
and to infer the properties of the early galaxies produc- 
ing it. In addition, its predictive power allows to design 
specific tests that could validate the results. 

The plan of the paper is the following. In §2 we present 
the N-body numerical simulations and in §3 we describe 
our assumptions about star formation. Then we turn to 
the discussion of radiative (§4) and stellar (§5) feedbacks 
and a brief sketch summarizing the possible galactic evo- 
lutionary tracks due to such processes is presented in §6. 
Results are given in §7 and further discussed in §9. A short 
summary concludes the paper. 

2. NUMERICAL SIMULATIONS AND HALO SELECTION 

We have simulated structure formation within a peri- 
odic cube of comoving length L = 2.55h~^ Mpc for a crit- 
ical density cold dark matter model {Clo=l, h=0.5 with 
a8=0.6 at z=0). Our choice for the normalization (cts) 
corresponds roughly to those inferred from the present- 
day cluster abundance (see, e.g.. Eke, Cole & Frenk 1996 
and Governato et al. 1999). We remind the reader that 
once the present-day value of o"8 for the SCDM run is se- 
lected, the redshift epoch of all other outputs with lower 
values of ag is uniquely specified. For example, for the 
case where the present-day normalization is chosen to be 
(78=0.6, CT8=0-3 output corresponds to the z = 1 epoch. It 
is useful to note that the rms linear overdcnsity of a sphere 
with a mass equal to the box is 2.1/(l-|-z). At early times 
when the amplitude of fluctuations are small for scales 
larger than the box size the simulation box should be large 
enough to be representative of the world model as a whole. 

^http;//www-hpcc.astro.washington.edu/tools/FOF/ 



But as the scale of non- linear fluctuations grows there will 
come a time when the simulation volume will not be large 
enough to provide a proper average. By redshift 8, for ex- 
ample, the rms linear overdensity of a sphere with the mass 
of the entire simulation is 0.23 and by then the mass spec- 
trum of dark halos in the simulation box may not match 
the real mass spectrum very well - particularly at the high 
mass end. 

The initial conditions for the simulations, were set up at 

z=100 and integrated in time to z=8.3 (with intermediate 
outputs at z=29.6, 25.3, 22.1, 19.8, 18.0, 16.5, 15.4, 14.3, 
11.4 and 10.9). We used a transfer function for CDM cal- 
culated with CMBFAST and a baryonic fraction i7b=0.06 
(Copi, Schramm & Turner 1995). The simulation, which 
uses 256'' particles (about 17 million), was computed with 
AP3M (Pearce & Couchman 1997), a parallel P3M code. 
A cubic spline force softening of 0.05 kpc Plummer 
equivalent was used so that the overall structure of the 
halos could be resolved. The good mass (particle mass 
is « 5 X 10'^ il/,;,) and force resolution allows us to study 
in detail the evolution of all structures likely to host the 
formation of primordial stars. 

In numerical simulations, halos can be identified using 
a variety of schemes. Of these, we have chosen one that is 
available in the public domain: FOF^ (Davis et al. 1985). 
In this scheme, all particle pairs separated by less than 
b times the mean interparticle separation are linked to- 
gether. Sets of mutually linked particles form groups that 
are then identified as dark matter halos. Other halo finders 
that are often used in literature to find virialized halos are 
HOP (Eisenstein & Hut 1998), the "spherical overdensity 
algorithm " or SO, that finds spherically averaged halos 
above a given overdensity (Lacey & Cole 1994) and the 
scheme recently developed by Gross et al. (1998). In the 
present study, we adopted the linking length that Lacey 
& Cole (1994) deduced to identify virialized halos with 
mean densities of «200 times the critical density at the 
epoch under consideration. Halo masses (and ultimately, 
our findings) do not depend on the halo finder used (EOF, 
HOP or SO), as long as it selects halos including all parti- 
cles down to the same overdensity. Systematic deviations 
between different algorithms are of the order of a few per 
cent (see Governato et al. 1999, Lacey & Cole 1994). In 
our analysis, we only consider halos consisting of 16 par- 
ticles or more, corresponding to a mass of 8 x 10^ Mq. 

We have compared the results of the numerical simula- 
tions for the dark matter halo distribution, with the one 
obtained from the Press & Schechter (1974, PS) formal- 
ism, so that we can use this semi-analytical expression 
to calculate certain quantities, such as the soft-UV back- 
ground (see § 4.2). The power spectrum normalization is 
a8=0.6, consistent with the above numerical simulations. 
Although PS slightly overestimates the halo number den- 
sity, the match is remarkably good over the entire range 
of redshifts considered. 

Finally, the correlation function for galaxies at high red- 
shift has been recently predicted (Governato et al. 1998) 
and then measured by Giavalisco et al. (1998), giving 
strong support to hierarchical models of galaxy formation. 
Bright galaxies at high redshift form on top of the highest 
peaks in the mass distribution, and so they are strongly 
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biased with respect to the underlying dark matter distribu- 
tion. In hierarchical models we expect a similar behavior 
for the most massive halos hosting Pop III stars. This mea- 
surement would give a strong constraint on the shape of 
the power spectrum at scales smaller than 1 Mpc, very rel- 
evant for all theories of galaxy formation. Indeed numcri- 
(;al simulations have already shown a large discrepancy be- 
tween the shape of the rotation curves of CDM halos and 
those of real dark matter dominated galaxies (see Moore 
et al. 1999a) as CDM models predict dark matter halos 
with steeper density profiles than those of real galaxies. 
We have preliminarly obtained the two-point correlation 
function of halos significantly contributing to the photon 
production between redshift w20 and 8.3. In our SCDM 
model the correlation hmction of halos hosting Pop III ob- 
jects has a (comoving) length scale of about 0.15 Mpc h~^, 
slowly decreasing with redshift. This decrease is expected 
(Tegmark & Peebles 1998) as the number of halos of mass 
sufficient to host Pop III objects increases with time. Feed- 
back effects of the; type discusscxl here are clearly strongly 
influencing the evolution of the correlation function. We 
plan to present the results concerning the correlation prop- 
erties of these halos in a future communication. 

3. FORMING THE FIRST STARS 

Once the gas, driven by gravitational instabilities, has 
been virialized in the potential well of the parent dark 
matter halo, further fragmentation of the gas and ignition 
of star formation is possible only if the gas can efficiently 
cool and lose pressure support. For a plasma of primordial 
composition at temperature T < 10'' K, the typical virial 
temperature of the early bound structures, molecular hy- 
drogen is the only efficient coolant. Thus, a minimum H2 
fraction is required for a gas cloud to be able to cool in 
a Hubble time. As the intergalactic relic H2 abundance 
falls short of at least two orders of magnitude with respect 
to the above value, the fate of a virialized lump depends 
crucially on its ability to rapidly increase its II2 content 
during the collapse phase. Tegmark et al. (1997) have ad- 
dressed this question in great detail by calculating the evo- 
lution of the H2 abundance for different halo masses and 
initial conditions for a standard CDM cosmology. They 
conclude that if the prevailing conditions are such that a 
molecular hydrogen fraction of order of /j^^ w 5 x 10"'' is 
produced, then the lump will cool, fragment and eventu- 
ally form stars. This criterion is met only by larger halos 
implying that for each virialization redshift there will ex- 
ist some critical mass, Merit , such that protogalaxies with 
total mass Al > Merit will be able to form stars and those 
with M < Merit will fail (see their Fig. 6 for the evolution 
of McriL with the virialization redshift). In reality, even 
halos with masses smaller than Merit could eventually col- 
lapse at a later time (Haiman & Loeb 1997), with a delay 
increasing with decreasing baryonic mass for a fixed rms 
amplitude of the fluctuation. However, as we will see be- 
low, these structures will be strongly affected by various 
feedbacks which essentially erase their contribution to the 
reionization process; hence we neglect this effect in our cal- 
culations. In the absence of additional effects that coifld 
prevent or delay the collapse (such as the feedbacks dis- 
cussed below), we can associate to each dark matter halo 
with M > Merit, a corresponding baryonic collapsed mass 
equal to Mf, — Q^M. Throughout the paper we adopt a 



value of the baryon density parameter Qi, — 0.06 consis- 
tently with the numerical simulations; note that the set 
of cosmological parameters used here are the same as in 
Tegmark et al. (1997), thus allowing a direct use of their 
results. 

As the gas collapses and the first stars form, stellar pho- 
tons with energies in the Lyman- Werner (LW) band and 
above the Lyman limit, respectively, can photodissociate 
H2 molecules and ionize H and He atoms in the surround- 
ing IGM. By this process a photodissociated/ionized re- 
gion will be produced around the collapsed object, whose 
size will depend on the source emission properties. The ra- 
diation spectrum, j(!^), adopted here is obtained from the 
recently revised version of the Bruzual & Chariot (1993, 
BC) spectrophotometric code, for a Salpeter initial mass 
function (IMF), a single burst mode of star formation, and 
a metallicity Z = 10~^Zq. This choice is supported by 
recent observational results concluding that the IMF in 
nearby systems has a slope close the Salpeter value above 
a solar mass, while fiattening at lower masses (see Scalo 
1998). Usually, a time independent IMF is assumed by 
most studies, but indirect evidences for an IMF biased to- 
wards massive stars at early times emerge (sec Larson 1998 
and references therein). Larson (1998) suggests that the 
IMF might be well represented by a universal power-law at 
large masses, flattening below a characteristic mass whose 
value is changing with time. The rationale for this con- 
clusion is that the mass scale for star formation probably 
depends strongly on temperature and since star-forming 
clouds were probably hotter at earlier cosmic times, the 
mass scale should also have been correspondingly higher. 
This results in an increase of the relative number of high- 
mass stars formed, and consequently of the number of ion- 
izing photons. The relevance of the Jeans mass scale for 
the IMF is still subject of lively debate and it is not clear 
to which extent it might be responsible for the mass dis- 
tribution. In view of these uncertainties and for sake of 
simplicity, we will adopt the Salpeter IMF throughout the 
paper. 

Fig. 1 shows the spectral energy distribution of an ob- 
ject per solar mass of stars formed at four different burst 
evolutionary times t = (0,1,2,32) x 10^ yr. At earlier 
times the spectrum is rather flat between the 13.6 eV and 
the Hell edge, with a pronounced cutoff at higher energies; 
at later stages, it becomes softer with strong discontinu- 
ities both at 13.6 eV and 24.6 cV. The number of ionizing 
photons decreases considerably as the stars age. 

The total luminosity of an object is determined 
by the mass of stars formed. For a baryonic mass 
M(,=10^ Mf,,5 Mq, the corresponding stellar mass is: 

M, - Mbfbf, ~ 1.2 X 10^ Mfa,5/b,8/*,i5 Mq, (1) 

where /(,=0.08 fb,8 is the fraction of virialized baryons that 
is able to cool and become available to form stars and 
fi,=0.15 /*.i5 is the star formation efficiency. With these 
assumptions the luminosity per unit frequency at the Ly- 
man limit, jo, as obtained from the adopted spectrum at 
early evolutionary times (Fig. 1), can be explicitly written 
in terms of Mf, : 

jo = 4 X 10^°M^ crgs^^Hz-^ 

~ 4.8 X 1023Mfe,5/MA,i5 ergs-'Rz-\ (2) 
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of which only a fraction jo/esc— 9-6x 10 M6^5/5.8/*.i5/e.sc,20 
erg s^^ Hz^^ (whore fpg,, = 0.2 /esc, 20 is the photon escape 
fraction from the proto-galaxy) is able to escape into the 
IGM. This accounts at least approximately for absorption 
occurring in the host galaxy. 

As fb, /* and fesc are the main parameters involved in 
the calculation, aside from the cosmological ones fixed by 
the simulations, it is useful to discuss them in more de- 
tail. Primordial stars can only form in the fraction of gas 
that has been able to cool, and the cooled-to-total bary- 
onic mass ratio in the collapsed objects is only a few per- 
cent. Numerical simulations (Abel et al. 1997b; Bromm, 
Coppi & Larson 1999) have shown that can be f«8%, 
which we assume as our fiducial value. The star formation 
efficiency, /^,, is rather uncertain and dependent over the 
properties of the star formation environment. Even its def- 
inition is not completely unambiguous, as star formation 
usually takes place in the densest cores of molecular clouds 
and the efficiency in the cores is obviously different from 
the one averaged on the parent molecular cloud. We define 

as the fraction of cooled (probably molecular) gas that 
has been converted into stars. The formation of a bound 
cluster system requires a star formation efficiency of nearly 
50%, where the cloud disruption is sudden and nearly 20% 
where cloud disruption takes place on a longer timescale 
(Marguhs & Lada 1983; Mathieu 1983). The highest ef- 
ficiency estimated for a star-forming region is about 30% 
for the Ophiuchi dark cloud (Wilking & Lada 1983; Lada 
& Wilking 1984). Lada, Evans II & Falgarone (1997) have 
studied the physical properties of molecular cloud cores 
in L1630, finding efficiencies ranging from 4% to 30% for 
different cores. Pandey, Paliwal & Mahra (1990) have in- 
vestigated the effect of variations in the IMF on the star 
formation efficiency in clouds of various masses. They con- 
clude that the efficiency is lower if the most massive, and 
consequently most destructive, stars form earlier. On a 
larger scale, Planesas, Colina & Perez-Olea (1997) have 
studied the molec;ular gas properties and star formation in 
nearby nuclear starburst galaxies, showing the existence 
of giant molecular clouds {Mh-2 ~ 10^ — IQ^Mq) with as- 
sociated HII regions where the star formation process is 
characterized by being short lived (< 3 x 10^ yr) and with 
an overall gas to stars conversion less than 10% of the gas 
mass. Giannakopoulou-Creighton, Fich & Wilson (1999) 
have recently determined the star formation efficiency in 
two MlOl giant molecular clouds, finding the values 6% 
and >11 %. In spite of great efforts to derive such an 
important quantity, various authors obtain results uncom- 
fortably different, even when studying the same star form- 
ing complex (see the prototypical case of 30 Dor). Given 
this situation, we use the educated guess f^, = 15%. Fi- 
nally, the value fesc=0-'2 is an upper limit derived from 
observational (Leitherer et al. 1995; Hurwitz, Jelinsky & 
Dixon 1997) and theoretical (Dove & Shull 1994; Dove, 
ShuU & Ferrara 1999) studies. The latter authors, in par- 
ticular, have included the effects of superbubbles in the 
computation of fesc, finding a value fesc ~ 6% for the 
case of coeval star formation history, which should be rel- 
evant to the present case. These studies concentrate on 
large disk galaxies like the Milky Way and it is not clear 
to which extent they can be applied to the small objects 
populating the early universe; also, the much lower dust- 
to-gas ratio could make the escaping probability higher. 



Again, fesc = 0.2 can only be seen as a reference value. 
It is worth noting that while for all processes considered 
here only the product fb x enters the model, fesc is in- 
stead an independent parameter; this reduces the number 
of effective free parameters to two. 

The total ionizing photon rate, 5^(0), from the formed 
stellar cluster can be written as: 

where i'/ = 13.6 eV and t/„=150 eV. As seen from Fig. 1 
the stellar spectrum drops off sharply above 100 eV even 
shortly after the burst; the choice of the upper integration 
limit follows from this remark. For the spectrum in Fig. 1 

at t = 0, eq. (3) becomes: 

SiiO) = 5.69 X 10~\jo/ hp) fesc 

~ 8.18 X 10^8 Mb,5fb,8f.,wfesc,20 8"^ (4) 

analogously we can define the LW photon production rate, 
Slw, which is typically Slw = PSi{0), with /3 « 1.5 at 
early evolutionary times. The UV photons create a cos- 
mological HII region in the surrounding IGM, whose ra- 
dius, Ri, can be approximated by the Stromgren (proper) 
radius, Rs = [3Si{0)/{ATmjja^^^)]^/^, where = 8 x 
10~^r2b/i^(l -|- z)^ cm~^ is the IGM hydrogen number den- 
sity and a'^^ is the hydrogen recombination rate to levels 
> 2. In general, Rg represents an upper limit for Ri, since 
the ionization front fills the time- varying Stromgren radius 
only at very high redshift, z « 100 (Shapiro & Giroux 
1987). For our reference parameters it is: 

Ri<Rs= 0.05 (Ofe/i') "^/^ (1 + z)^^ Sl(^ kpc, (5) 

where 847 = S^{0)/{10'^'^ s-^) and (1 + 0)30 = (1 + z)/SO. 
In our calculation we use the exact solution for Ri as nu- 
merically calculated in CFA, and whose analytical approx- 
imation is given by Shapiro & Giroux (1987). Typically, 
Ri is about 1.5-2.0 times smaller than Rs and cosmologi- 
cal expansion cannot be neglected. Clearly, these solutions 
assume that the ionizing object is irradiating a volume of 
the IGM which was previously neutral. If instead the ob- 
ject is located inside (or overlapping) a preexisting ionized 
sphere, a detailed solution of the radiative transfer would 
be necessary. This occurrence can result in an underes- 
timate of our calculated ionized volume when substantial 
sphere overlapping takes place, i.e. close to the reioniza- 
tion epoch. 

Similarly to Rs, we can define the Hell Stromgren ra- 
dius as Rs^He — [3S'i,_f/e(0)/(47rn^e++neQ://e++)]"'^^^j where 
Si,He{0) is the rate of ionizing photons with hiy > 54.4 eV, 
nHe++ and rig are the He^^ and electron number density, 
respectively, and aHe++ is the He++ recombination rate. 
For the adopted spectrum it is: 

SiM^) = 1-87 X 10-Hjo/hp)fesc 

~ 2.68 X 10^6 Mb,5fb,sf*,15fesc,20 s'^ (6) 

We have derived nHe++ and Ue assuming that the helium 
and hydrogen neutral fraction is negligible. For our refer- 
ence parameters and given eqs. (4) and (6), the radius of 
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the sphere of doubly ionized hehum is: 



Ri,He ^ Rs,He = 6 X 10"^ (l^feft^) (1 + 2)"^ 



ol/3 
30 "^47 



kpc, 
(7) 

thus RsMe ^ O.IR,. 

Given a point source that radiates S lw photons per sec- 
ond in the LW bands, a good estimate of the maximum 
radius of the H2 photodissociated sphere is the distance 
at which the (optically thin) photo-dissociation time be- 
comes longer than the Hubble time: 



Rd ^ 2.5 /i- 1/2(1 + z)^^'^S]!^^^, kpc, 



(8) 



where S'lw,47 = 5'lm'/(10**^ s"^). CFA also considered the 
time dependent evolution of both the ionization and disso- 
ciation fronts; here we use the exact value for i?j and the 
equilibrium values given by the equations above for both 
Ri^He and R^- The latter values represent upper limits 
to the extent of the spheres as emphasized above: ionized 
regions at high rcdshift fill only partially the correspond- 
ing Stromgren spheres, while for the dissociated regions 
the radius Rd is completely filled, although on a relatively 
long time scale. For example an object born at z = 30 will 
have a fully developed surrounding dissociated sphere by 
z « 20. 

To each dark matter halo in which stars can form we 
assign the ionized and dissociated spheres produced by its 
stellar cluster according to the above prescriptions. This 
allows to derive the three dimensional structure and topol- 
ogy of the ionized/dissociated regions as a function of cos- 
mic time. 

4. RADIATIVE FEEDBACK 

In addition to their local effects, the first objects will 
also produce UV radiation which could in principle in- 
troduce long range feedbacks on nearby collapsing halos. 
Particularly relevant is the soft UV background in the LW 
bands, as by dissociating the H2, it could influence the star 
formation history of other small objects preventing their 
cooling. 

CFA have argued that the SUVB produced by Pop Ills 
is below the threshold required for the negative feedback 
on the subsequent galaxy formation to be effective before 
z « 20. At later times this feedback becomes important 
for objects with T„jr ~ Th = 10000 K, corresponding to a 
mass Mh=4:A x lO^M© (1 + Zy.,ry^-^h^^ (Padmanabhan 
1993), where Zyir is the redshift of virialization, for which 
cooling via Lya line radiation is possible. Proto-galaxies 
with M > Mff are not affected by the negative feedback 
and are assumed to form. In principle even for these larger 
objects a different type of feedback can be at work: halos 
with a circular velocity lower than the sound speed of the 
gas in an ionized region may be prevented from collapsing 
due to gas pressure support in the gravitational potential. 
Thoul & Weinberg (1996) (see also Babul & Rees 1992) 
have shown that the collapse is only delayed by this pro- 
cess, and therefore we neglect this complication. 

At higher redshift the radiative feedback can be induced 
by the direct dissociating flux from a nearby object. In 
practice, two different situations can occur: i) the col- 
lapsing object is outside the dissociated spheres produced 
by preexistent objects: then its formation could be af- 
fected only by the SUVB {JLW,b), as by construction the 



direct flux (Jlwm) can only dissociate molecular hydro- 
gen on time scales shorter than the Hubble time inside 
Rd; ii) the collapsing object is located inside the disso- 
ciation sphere of a previously collapsed object: the ac- 
tual dissociating flux in this case is essentially given by 
JLW.max = iJLW,b + JLW,d)- It is thus assumed that, given 
a forming Pop III, if the incident dissociating flux {JLW,b 
in the former case, JLW,max in the latter) is higher than 
the minimum flux required for negative feedback (Jg), the 
collapse of the object is halted. This implies the existence 
of a population of "dark objects" which were not able to 
produce stars and, hence, light. 

4.1. Minimum flux for negative feedback 

As already stated, in the absence of an external dissoci- 
ating fltix, an object can form if it satisfies the condition 
M > Merit- In this subsection we investigate how this 
minimum mass is changed by the requirement that the 
object is able to self-shield from an external incident dis- 
sociating flux that could deplete the H2 abundance and 
prevent the collapse of the gas. To this aim, we consider 
the non-equilibrium multifrequency radiative transfer of 
an incident spectrum inside a homogeneous gas layer, and 
study the evolution of the following nine species: H, H^, 
H+, He, He+, He++, H2, Hj and free electrons. The 
energy equation can be written as (Ferrara & Giallongo 
1996): 



k 



(7 - 1) dt 



X 



nx 



where T, n and p are the gas temperature, number density 
and pressure, respectively. The gas is assumed of primor- 
dial composition with a helium abundance equal to nne = 
0.1 uh and total number density n. The symbol X= H, 
He, H2 denotes the species considered, and i its state of 
ionization; obviously, i = for hydrogen, and i = 0, 1 
for helium; hx = nx/n is the relative abundance of the 

species X\ fix- = {i-max + 1)^[1 + X]-""'"(« + l)xxi+i]~\ 
with imax = 0, 1 for H and He, respectively and Xxi+^ = 
nxi+i /nx the fractional density of the ionization state i+1 
of the element X; 7 is the specific heat ratio assumed to 
be constant and equal to 5/3. The function C represents 
the net cooling rate per unit volume and is given by the 
difference between heating and cooling: 



^X') = Y.I1 <X')H{X') - 5] 5^ A(X^ 



(10) 



X 



The cooling term A (in erg cm~'^ s"^) includes colli- 
sional ionization and excitation of H, He, He"*", recombina- 
tion to H, He, He+, dielectronic recombination to He, free- 
free (Black 1981), Compton cooling (Peebles 1971) and H2 
cooling (Martin, Schwarz & Mandy 1996); Ti includes the 
heating terms due to photoionization of H, He, He+, H2, 
photodissociation of H2' and H2, and H~ photodetachment 
and is given by: 



' 4TTJs{hiy) 

dv h{v - Vx)(Jxi{v), (11) 



where vx indicates the ionization limit for each species, a 
is the photoionization cross-section and Jg is the incident 
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flux, given below. The various cross sections are given in 
Abel et al. (1997b), apart from (J21 (taken from Brown 
1971) and 0-25 (Shapiro & Kang 1987), and they are num- 
bered according to the nomenclature of Abel et al. (1997b). 
We have assumed that He, and are in equilibrium 
as all reactions determining the abundance occur on 
much shorter time scales than those relevant to the H2 
chemistry; He and H2' do not influence substantially the 
final results. The chemical network includes the 27 re- 
actions listed in Abel et al. (1997b). We have adopted 
the same rates of that paper except for ki (taken from 
Cen 1992), k2 (Black 1981), ku and k^ (Shapiro & Kang 
1987). The photoionization rate 7j, is given by: 

7,(X0 = 1^°° di^^^^^^aiX^) [1 + HXl] , (12) 

where cf) is the secondary ionization rate. We have adopted 
the "on the spot" approximation in which the diffuse field 
photons are supposed to be absorbed close to the point 
where they have been generated. 

To derive generally valid results for the negative feed- 
back we have taken the same incident spectrum presented 
in Fig. 1 but at a later evolutionary time (2 x 10'' yr after 
the burst) in order to roughly average over the different 
population evolutionary stages. An useful analytical ap- 
proximation to the adopted spectrum is: 

{f3 E < 13.6 eV, 

{E/n.&y^-^ 13.6 <£;<25eV, 

5 X io-3e-(i=-25)/3.3 25 <E<m eV, 
2.6 X io-6e-(^-50)/i-5 50 < s < 15O eV, 

(13) 

with /3=50 and E = hv. Here Jg^ (units of erg s ^ cm ^ 
Hz~-^ sr~^) is the parameter with respect to which we 
quantify the negative feedback at a given redshift. 

Initially, we assume that the gas layer is homogeneous, 
with a mean density (pt) 18 tt^ times higher than the 
background intergalactic medium. The fraction of free 
electrons is taken to have the relic value after recombi- 
nation X = lo-^vi^^ny^h''^ ( see for example Blumenthal 
et al. 1984); in principle, x could have changed during the 
virialization process: we neglect this possibility. Wc use 
equilibrium values for the He ionic abundances; H~ and 
H2' abundances can be arbitrarily small, as their choice 
does not influence the results. The initial condition for 
the molecular abundance is slightly more delicate; we de- 
termine it as follows. For a given Merit (the mass of an 
object able to collapse in the absence of radiation) we can 
determine the corresponding virial temperature Tyir . This 
represents the initial condition for the gas temperature. 
The initial H2 abundance is set by the value required to 
satisfy the condition that the free-fall time of the object is 
equal to its cooling time at temperature T^r ■ This assures 
that, if no feedback is imposed on the object, it would col- 
lapse and form stars. 

Starting from these initial conditions we have followed 
the chemical evolution up to a free-fall time. These condi- 
tions are such that in the absence of radiation the fraction 
of H2 is high enough to allow for the collapse of the ob- 
ject on a time scale comparable to the free-fall time. The 
effect of the radiation field is to decrease this abundance 
in the external regions which therefore cannot cool. In 



the interior, instead, where the H2 abundance essentially 
remains equal to the initial one, the collapse can proceed 
unimpeded. For this reason the collapse conditions must 
be checked after a free-fall time. As an example we dis- 
cuss the case of an object which virializes at z = 29.6 and 
to which an incident flux with spectrum given by eq. (13) 
with Jsfi = J2ilO~^^ erg s~^ cm~^ Hz~^ sr^-^ has been im- 
posed. The critical mass for collapse. Merit-, corresponding 
to that redshift is Merit = 9 x 10^ Mq ; this fixes the initial 
H2 fraction as explained above. Figs. 2-4 show the cor- 
responding hydrogen ionization fraction, x, the gas tem- 
perature, T, and the neutral molecular hydrogen fraction, 
/jj2 5 after a free-fall time. For sake of comparison, we also 
show the cases relative to a power law spectrum (PL) with 
q;=1.5 and a cutoff energy of 40 keV. Different values of the 
(3 parameter, describing the relative number of photodisso- 
ciating and ionizing photons in the spectrum, are studied 
in order to cover a wide range of evolutionary phases of 
the stellar energy emission. As seen from Fig. 1, the value 
of (3 tends to increase as the stellar chistcr ages. Our stan- 
dard case is the BC spectrum with /3=50, appropriate to 
an age of 2 x 10^ yr. 

From Fig. 2, it is clear that the ionization fraction is 
mostly determined by photons with energies above the 
Lyman limit; however, in the cases with /3=1, as the fi- 
nal abundance of H2 is increased with respect to the ini- 
tial one (see Fig. 4), electrons are depleted by the H~ 
molecular hydrogen formation channel, with respect to 
the cases with /3=50. The high energy photons present 
in the power-law spectrum penetrate deeper inside the 
layer and consequently the initial ionization is enhanced 
over a larger depth than for the BC spectrum, although 
with a lower ionization level; this feature is typical of hard 
ionizing spectra. For the same reason, a PL spectrum 
heats the cloud to higher temperatures and the plateau 
extends into deeper regions (Fig. 3). Lower temperatures 
and higher fn-i are obtained for smaller values of /3 (see 
Fig. 3), when the H2 cooling is more efficient. It is clear 
from Fig. 3, that, although in the external regions of the 
layer heating is dominant due to photoelectric effect, in the 
inner regions, where the medium becomes optically thick 
to LW photons, after a free-fall time the temperature has 
decreased by about 30% due to H2 cooling. In Fig. 4 the 
profile is shown. In the external regions, the most 
important mechanisms for the H2 destruction are direct 
photoionization by photons with energies above 15.4 eV 
and LW photons dissociation. A PL spectrum with the 
same value of (3 yields a larger H2 destruction, while if 
the upper PL spectrum cutoff is increased the curves shift 
towards larger depths due to penetrating photons. Note 
that when /3=1, H2 is formed rather than destroyed due to 
the larger abundance of free electrons and paucity of LW 
photons available. Finally, the fH^. smoothly approaches 
the initial value in the internal regions due to the loss of 
free electrons (see Fig. 2). 

We would like to briefly comment on the differences be- 
tween the above curves for the H2 abundance and the anal- 
ogous ones derived by HRL. These are mainly due to: (i) 
a different incident flux and (ii) their assumption of chem- 
ical and ionization equilibrium at a constant temperature, 
which could not be reached within a free-fall time. As they 
also note in a previous paper (Haiman, Rees & Loeb 1996), 
this might lead to an overestimate of the H2 abundance. In 



8 



addition, they used the Lepp & ShuU (1984) eoohng func- 
tion, that, at the low densities of interest here, is higher 
than the one we have adopted (for a comparison see Galli 
& Palla 1998). 

As shown above, after a free-fall time the temperature 
in the inner regions of the protogalactic cloud has con- 
siderably decreased due to H2 cooling. It is normally as- 
sumed that this process would eventually lead to the col- 
lapse, fragmentation and star formation. These regions are 
therefore not affected by the imposed external flux as they 
can efficiently self-shield. The gas in the outer regions, 
up to the shielding radius, Rsh{Js.o), is instead globally 
heated by radiation and fails to collapse. On this basis, 
we can then define the minimum total mass required for 
an object to self-shield from an external fiux of intensity 
Jsfl at the Lyman limit, as Msh = (4/3) 7r(p/j)i?^^, where 
{ph) = {pb)fl'^^ is the mean dark halo matter density. In 
practice, Rsh is defined as the point where the temper- 
ature profile has increased by 0.01% from the inner flat 
curve (see Fig. 3). Values of Mgh for different values of Jgfi 
have been obtained at various redshifts and will be then 
used in the calculations presented below; these curves are 
shown in Fig. 5. We point out that the results depend only 
very weakly on the detailed shape of the spectrum as long 
as this is produced by stellar sources with a relatively soft 
spectrum; differences up to a factor of several can be found 
if instead a hard spectral component is present (Haiman, 
Abel & Rees 1999). In this paper we restrict our analysis 
to the first case. Protogalaxies with masses above Mh for 
which cooling is predominantly contributed by Lya line 
cooling will not be affected by the negative feedback stud- 
ied here: these objects lie on the upper dashed portions 
of the curves in Fig. 5. The collapse of very small objects 
with mass < Merit is on the other hand made impossible 
by the cooling time being longer than the Hubble time 
(lower dashed portion of the curves). Thus the only mass 
range in which negative feedback is important (solid por- 
tion) lies approximately in 10^ — IO^Mq, depending on 
redshift. In order for the negative feedback to be effective, 
fluxes of the order of 10"^^ - lO'^^erg s"^ cm'^ Hz"^ 
sr^^ are required. In turn, by using eq. (4), we see that 
these fluxes are produced by a Pop III with baryonic mass 

1 /2 

M6^5 at distances closer than ~ 21 — 7 x Mj, 5 kpc for the 
two above flux values, respectively, while the SUVB can 
reach an intensity in the above range only after z ^ 15 
(see § 7.2). This suggests that at high z negative feedback 
is driven primarily by the direct irradiation from neigh- 
bor objects in regions of intense clustering, while only for 
z ^15 the SUVB becomes dominant. 

4.2. Maximum incident flux 

We define the maximum dissociating flux felt by a form- 
ing Pop III as the sum JLW,max = {JhWfi + JLW,d), where 
Jhwfi is the SUVB and JLW,d is the direct flux produced 
by a nearby luminous object. JLW,d is defined as: 



'LW,d = 



475-7-2 ' 



(14) 



r is the distance between the forming object and the source 
of the dissociating fiux. The value of jo is the one given in 
eq. (2) and the value of /3 = 1.5 is adopted for the reasons 
explained in §3. 



We now derive the intensity of the SUVB. As the LW 
range is very narrow (w 2 eV) we consider the average flux 
at the central frequency of the band, hvo = 12.45 eV. To 
properly calculate the intensity of the SUVB we must con- 
sider the intergalactic H2 attenuation including the effects 
of cosmological expansion and treat in detail the radiative 
transfer through LW II2 lines. LW lines are optically thin, 
with Ti = NH2,i<^i ~ 0.05, essentially for all the lines. This 
implies that as Jlw is redshifted due to cosmological ex- 
pansion through LW lines, it is attenuated by each line by 
a factor e""^'. Abgrall & Roueff (1989) have included in 
their study of classic H2 PDRs more than 1000 LW lines. 
As in our case H2 formation, which leaves the molecule 
in excited roto/vibrational levels, is negligible, a smaller 
number of lines w 70 involving the ground state only - 
needs to be considered. Globally, Jlw is attenuated by a 
factor e""^"^, where th2 = "Y^iTi, and i runs up to the 71 
lines considered from the groimd roto/vibrational states; 
we obtain th^ ~ 3, depending on the number of lines en- 
countered, and thus on the photon energy. The intensity 
of the background is: 



JLW,b{z) = c [ e{iy\ z') Ah^h^ TTT"^ 
Jz„„. 0- + Z' 



dt 
dz' 



e{v',z') 



JMmi 



3{v')N{M, z')dM, 



dz', (15) 



(16) 



where z,. 



29.6 corresponds to the redshift at which 



the first objects start to form, u' = vo{l + z')/{l + z); 
e{u' , z') [erg cm~^ s~^ Hz~^] is the proper emissivity due 
to sources with M^in < M < Mmax; Ah.h-i takes into 
account the H and II2 line absorption in the LW band and 
the H absorption above the Lyman limit; the correspond- 
ing curves of growth are taken from Federman, Glassgold 
& Kwan (1979). is taken from the BC spectrum 

as described in § 3 with intensity given by eq. (2) and 
Af{M, z') is the source number density in the mass inter- 
val dAd derived from the PS formalism. The upper cut-off 
is Mmax = IO^^Mq, the maximum statistically significant 
mass down to z =8.3, the lowest redshift of the numer- 
ical simulation. Mmin is the mass of the smallest halo 
that has been able to collapse at the previous redshift 
step; we are thus assuming that all objects with masses 
greater than Mmin will contribute to the SUVB. This cal- 
culation of Jlwm is correct as long as the filling factor of 
the ionized and dissociated regions is small. As the vol- 
ume filling factor of these regions grows with time, their 
attenuation of the radiation emitted by the sources is di- 
minished and eventually, when the dissociation/ionization 
process is complete, vanishes. To roughly take into ac- 
count these complications, we assume that only a fraction 
1 — /i of the IGM gas, where fi is the filling factor of either 
ionized or dissociated regions, contributes to the attenua- 
tion. We note that taking into account this effect, together 
with the determination of the lower limit Mmin of the in- 
tegral governed by the negative feedback, implies that eqs. 
(15)-(16) cannot be solved directly but must be computed 
iteratively as we derive the evolution of the reionization 
process. 

5. STELLAR FEEDBACK 
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Once the first stars have formed in the host protogalaxy, 
they can deeply influence the subsequent star formation 
process through the effects of mass and energy deposition 
due to winds and supernova explosions. While low mass 
objects may experience a blowaway, expelling their entire 
gas content into the IGM and quenching star formation, 
larger objects may instead be able to at least partially re- 
tain their baryons. However, even in this case the blowout 
induces a decrease of the star formation rate due to the 
global heating and loss of the galactic ISM. These two 
regimes are separated by a critical mass, M^j,, to be cal- 
culated in the following. 

For the relatively small objects present during the reion- 
ization epoch 30 ^ 2 <^ 10 the importance of these stellar 
feedbacks can hardly be overlooked. To understand the 
role of stellar feedback, let us consider a collapsed object 
with mass lower than Mby. Then the star formation is 
suddenly halted as the entire gas content is removed. In 
this case, the ionizing photon production will last only for 
a time interval of order toB ~ 10^ yr, the mean lifetime 
of the massive stars produced initially; after this, the ion- 
ized gas around the source will start to recombine at a 
fast rate as a result of the highly efflcient high z Compton 
cooling, rapidly decreasing the temperature inside the ion- 
ized region. Due to the short lifetime and recombination 
time scales, these object will only produce transient HII 
regions which will rapidly disappear. In fact, the recom- 
bination time scale, when the Compton cooling is taken 
into account, is of order tree ~ (1 — 50) x 10® yr at redshift 
2;Ri30-10, respectively, and therefore much shorter than the 
corresponding Hubble time. The dissociating photon pro- 
duction will nevertheless continue for a longer time, due to 
the important contribution of long-lived intermediate mass 
stars formed in the same initial star formation burst. This, 
combined with the fact that there is no efficient mechanism 
available to re-form the destroyed H2 in the IGM analo- 
gous to H recombination, implies that the dissociation is 
not impeded by blowaway and the contribution from these 
small objects should be yet accounted for. However, this 
is not necessary. In fact, after blowaway, H2 is efficiently 
formed in the shocked IGM gas, cooling under non equi- 
librium conditions (Ferrara 1998). The final radius of the 
cooled shell behind which H2 is formed by this process is: 



-19/10 



kpc. 



(17) 



We notice that the formation process produces an amount 
of H2 roughly similar to the one destroyed by radiation 
(Ferrara 1998). As inside Rs H2 is re-formed, the net ef- 
fect on the destruction of intergalactic H2 of blown away 
objects is negligible (if not positive) and we can safely ne- 
glect them in the subsequent calculations. 

In conclusion, low mass objects with mass < Mi,y pro- 
duce ionization regions which last only for a recombination 
time; larger objects can survive but their star formation 
ability is impaired by ISM loss/heating. The rest of the 
Section is devoted to quantify these two cases. 

5.1. Low m,ass objects: M < M^y 

Due to the small mass of these proto-galaxies, mas- 
sive star formation and, particularly, multi-supernova ex- 
plosions result in a blowaway. The blowaway is always 
preceded (unless the galaxy is perfectly spherical) by a 



blowout, as explained in FT (see also MF). The blowout 
and blowaway conditions can be derived by comparing the 
blowout velocity, vj,, and the escape velocity of the galaxy, 
Ve- Prom an analysis of the Kompaneets (1957) solution 
for the propagation of a shock produced by an explosion in 
a stratified medium, it follows that the shock velocity de- 
creases down to a minimum occurring at about 3H, where 
H is the characteristic scale height, before being reaccel- 
erated; Vb corresponds to such a minimum. 

To calculate these two velocities we define a proto- 
galaxy as a two component system made of a dark matter 
halo and a gaseous disk. We assume a modified isother- 
mal halo density profile ph{r) = Pe/[^ + (t'/^'o)^] extend- 
ing out to a radius rh = [3Mh/ {4:7r{18TT^)perit)]^^^ , de- 
fined as the characteristic radius within which the mean 
dark matter density is IStt^ times the critical density 
p^^it = 3F2Qq(i + zf/S-wG = 1.88 x IQ-'^^h'^ (1 + zf 
g cm~^; Mfi is the halo mass. For such a halo the escape 
velocity can be written as: 



Ve = (2|$(r,)|)V2 {AnpGpevl)'/' = {^-^ 
~ 0.9M^y/(l + zf '^]!iia s-\ 



1/2 



(18) 



where ^{vd) is the the gravitational potential of the halo 
(in this calculation we neglect the contribution of baryons 
to the potential) calculated at the radius of the disk, r^; 
p=1.65; the assumption ^ tq has been made. 

The explicit expression for Vb has been obtained by FT: 



Vb = 2.7 



L 



HHpd 



1/3 



km s 



(19) 



where L is the mechanical luminosity and {pd) is the disk 
mean gas density. L at redshift z can be written as: 



iff 



~ 9.4 X 10^5(1 +2)3/2M6,5/6,8/*,i5ergs s-1, 

where eo = 10° ergs is the energy of a SN explosion; as 
previously stated, we assume a Salpeter IMF, according to 
which one supernova is produced for each 135 Mq = 
of stars formed; the free-fall time is tff = {A-KG{ph))~^^'^\ 
Due to the dissipative nature of the gas, the baryons 
which should be initially distributed approximately as the 
dark matter, will lose pressure and collapse in the gravita- 
tional field of the dark matter. If the halo is rotating, the 
gas will collapse in a centrifugally supported disk and , as- 
suming that all the baryons bound to the dark matter halo 
are able to collapse in the disk, {pd) — {r\Ub/r'\H){ph) ■ 
The radius of the disk can be estimated by imposing that 
the specific angular momentum of the disk, j^, is equal to 
the halo one, jh (Mo, Mao & White 1998; Weil, Eke & 
Efstathiou 1998): 



jh = V^XVeVh = ^V, 



■ Jd, 



(21) 



where Vc is the circular velocity, id is the disk scale length 
and A the standard halo spin parameter. As shown by 
numerical simulations (Barnes & Efstathiou 1987; Stein- 
metz & Bartelmann 1995) A is only very weakly dependent 
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on Mh and on the density fluctuation spectrum; its dis- 
tribution is approximately log-normal and peaks around 
A = 0.04. From eq. (21) we then obtain ^c^ = {\/V2)rh. 
For an exponential disk, as implicitly assumed in deriving 
eq. (21) above, the optical radius (i.e. the radius encom- 
passing 83 % of the total integrated light) is 3.2 Id] then, 
the HI radius in galaxies is typically found to be « 2 times 
larger than the optical radius (Salpeter & Hoffman 1996). 
Thus, for the radius of the gaseous disk we take: 

Td = A.bXrh. (22) 

The scale height, H, is roughly given by: 

where Cg is the effective gas sound speed, including also a 
possible turbulent contribution; we take ~ 1 km s~^, 
corresponding to the minimum temperature allowed by 
molecular cooling. It is useful to note that H^{pd) = 
0.24:c'^Qb/G, i.e. is independent of mass and redshift. We 
point out that these estimates are in excellent agreement 
with the properties of Pop III disks found in numerical 
simulations of Bromm, Coppi & Larson (1999). 

In the following we derive the necessary condition for 
blowaway to occur. Following blowout, the pressure in- 
side the hot gas bubble drops suddenly due to the inward 
propagation of a rarefaction wave. The lateral walls of the 
shell, moving along the galaxy major axis will continue to 
expand unperturbed until they are overcome by the rar- 
efaction wave at rc, corresponding to a time tc elapsed 
from the blowout. After that moment, the shell enters the 
momentum-conserving phase, since the driving pressure 
has been dissipated by the blowout. The requirement for 
the blowaway to take place is then that the momentum 
of the shell (of mass Mc at Vc) is larger than the momen- 
tum necessary to accelerate the gas outside Tc at a velocity 
larger than the escape velocity: 

McVr^ > MhVe. (24) 

With the assumptions and algebra outlined in FT, one can 
write the condition for blowaway to occur: 

^ > (e-a)2a-V/2, (25) 

Ve 

where Vb is the blowout velocity, Ve is the escape velocity, 
a is a parameter equal to 2/3 and e = rd/H is the ratio 
of the major to the minor axis of the proto-galaxy. From 
eqs. (22) and (23), we find: 

e-' = 15^3A(fi)\ (26) 

Substituting eqs. (18), (19) and (26) into eq. (25), we find 
the critical mass below which blowaway occurs. It is use- 
ful to recall that this prescription has been carefully tested 
against hydrodynamical simulations presented in MF. 

In Fig. 6 Mby and Merit are plotted for the usual cos- 
mological parameters {ilo=l, h=0.5 and ^6=0. 06) and for 
the three runs A, B, C, whose parameters are defined in 
Tab. 1 and described in detail in §7. As Mby depends 



on the product fb x all other runs give the same re- 
sults as run A. As it is seen from the Figure, the blowaway 
mass depends very weakly on this product (approximately 
(X {fbf*y^'^), and it is always larger than Merit, but the 
number of objects suffering blowaway decreases with red- 
shift. Note as that the vast majority of the objects which 
virialized at redshifts XL 15 will be blown away due to their 
low mass predicted by hierarchical models. We also point 
out that at redshifts below w6, when objects with masses 
> IO^Mq start to dominate the mass function, our esti- 
mate of Mby becomes less accurate. This is due to the 
assumption implicitly made that all SNe drive the same 
supcrbubble. AlthoTigh valid for small objects, for larger 
galaxies this idealization becomes increasingly poorer and 
one should consider the fact that in general SNe are dis- 
tributed among OB associations with different values of L 
spread around the galactic disk. However, our calculations 
are not affected by this problem as we stop the simulations 
at z « 8. 

5.2. Larger objects: M > Mby 

Objects with M > Mby will not experience blowaway 
and rather than explode, they continue to produce stars 
quiescently. Nevertheless, a blowout with consequent mass 
loss and heating of the ISM by SN shock waves will very 
likely take place. MF found that galaxies with masses 
lower than a few x lO^^M© will suffer losses in an outflow, 
although the gas loss efficiency is not very high. Most 
important, perhaps, is the evaporation and destruction of 
molecular clouds as a result of an enhanced star formation 
activity which depresses the star formation rate. These 
processes are not yet fully understood in their complex- 
ity and detailed hydrodynamical models implementing the 
physics regulating the multi-phase structure of the ISM 
are required to make further progress. We parameterise 
these processes by an heuristic approach prescribing that 
the mass transformed into stars is lower than the value 
used for low mass objects by a mass dependent factor 

^feed' 

Mi = SfeedM^ = .,3 M^, (27) 

* ^ l + {Mby/M)t ' ^ ' 

where 7 is a free parameter usually taken equal to 2. We 
have checked that our results are almost insensitive to dif- 
ferent values of 7. This parameterization of the stellar 
feedback is analogous to the one typically used in semi- 
analytical models of galaxy formation (White & Frenk 
1991; Kauffman 1995; Baugh et al. 1998a; Guiderdoni et 
al. 1998). This factor allows for a maximum 50% decrease 
of the stellar mass of a galaxy due to gas heating/loss; 
this stellar feedback becomes increasingly less important 
for larger objects. 

6. SUMMARY OF EVOLUTIONARY TRACKS 

Given the large number of processes so far discussed 
and included in the model, it is probably worth to sum- 
marize them briefly before presenting the results. Fig. 7 
illustrates all possible evolutionary tracks and final fates 
of primordial objects, together with the mass scales de- 
termined by the various physical processes and feedbacks. 
We recall that there are four critical mass scales in the 
problem: (i) Merit, the minimum mass for an object to 
be able to cool in a Hubble time; (ii) Mh, the critical 
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mass for which hydrogen Lya Une coohng is dominant; 
(iii) Mgfi^ the characteristic mass above which the object is 
self-shielded, and (iv) M^y the characteristic mass for stel- 
lar feedback, below which blowaway can not be avoided. 
Starting from a virializcd dark matter halo, condition (i) 
produces the first branching, and objects failing to sat- 
isfy it will not collapse and form only a negligible amount 
of stars. In the following, we will refer to these objects 
as dark objects. Protogalaxies with masses in the range 
Merit < M < Mh are then subject to the effect of ra- 
diative feedback, which could either impede the collapse 
of those of them with mass M < Mgh, thus contributing 
to the class of dark objects, or allow the collapse of the 
remaining ones (M > Mgh) to join those with M > Mh 
in the class of luminous objects. This is the class of ob- 
jects that convert a considerable fraction of their baryons 
in stars. Stellar feedback causes the final bifurcation by 
inducing a blowaway of the baryons contained in luminous 
objects with mass M < M^y-, this separates the class in two 
subclasses, namely "normal" galaxies (although of masses 
comparable to present day dwarfs) that we dub gaseous 
galaxies and tiny stellar aggregates with negligible traces 
(if any) of gas to which consequently we will refer to as 
naked stellar clusters. The role of these distinct popula- 
tions for the reionization of the universe and their density 
evolution will be clarified by the following results. 

7. RESULTS 

In this Section we present the results obtained by in- 
cluding all the effects discussed above in the calculation of 
the reionization history. As the cosmological model (CDM 
with 0,0=1, h=0.5, flb=O.OQ) is fixed by the N-body simu- 
lations, there are three free parameters left, namely fi,, 
and /esc) defined in § 3 (actually, as already pointed out, 
the number of effective free parameters is two, /& x f^, and 
fesc)- We have performed runs with different values for 
these parameters, labelled as runs A-D in Table 1; our ref- 
erence case is run A (/;, = 0.08, = 0.15 and fesc = 0.2). 
We have made runs with the same parameters as in run 
A, but in which only objects with masses M > Mh are 
considered (run Al) and the normalization of the fluctua- 
tion spectrum has been varied (run A2). These cases are 
intended to isolate the effects of negative feedback and of 
a different evolution of the dark matter halos, as explained 
below. 

The present approach, based on a combination of very 
high resolution N-body numerical simulations and detailed 
treatment of the radiative/stellar feedbacks, allows us to 
study the inhomogeneous ionization process and the na- 
ture of the sources producing it at a remarkable detail 
level. The main output of our model consist of simulation 
boxes in which ionized/dissociated regions can be clearly 
identified and for which we can derive the global three 
dimensional structure as a function of cosmic time. 

7.1. Filling factor 

As a typical example, we show in Fig. 8 the actual out- 
come for run A. There, the spatial distribution of reion- 
izcd regions is shown at four different redshifts [z = 
22.1,19.8,18.0,15.4). The spheres grow in number and 
volume with time, although from z the number of 
luminous objects flattens. The reason for this behavior 
is that the population is dominated by small mass ob- 



jects whose collapse is prevented either by the condition 
M > Merit or by the increased intensity of the SUVB (see 
§ 7.4). The images shown in Fig. 8 give an immediate and 
qualitative view of the evolution of the volume occupied by 
ionized atomic hydrogen. A more quantitative measure of 
the filling factor of the ionized gas can easily be obtained. 

The filling factors of the dissociated II2 and ionized H 
are defined as the box volume fraction occupied by those 
species. The results are shown in Fig. 9a and 9b, for dif- 
ferent runs. The intergalactic relic molecular hydrogen is 
found to be completely dissociated at very high redshift 
{z « 25) independently of the parameters of the simula- 
tion. This descends from the fact that dissociation spheres 
are relatively large and overlap at early times. Ionization 
spheres are instead always smaller than dissociation ones 
and complete reionization occurs considerably later. Ex- 
cept for run C, when reionization occurs by primor- 
dial galaxies are able to reionize the IGM at a redshift 
z«10. The filling factor is approximately linearly propor- 
tional to the number of sources but has a stronger (cu- 
bic) dependence on the ionization sphere radius. Thus, 
although the number of (relatively small) luminous ob- 
jects present increases with redshift, the filling factor is 
only boosted by the appearance of larger, and therefore 
more luminous, objects which can ionize more efficicintly. 
The subsequent flattening at lower redshifts, is obviously 
due to the fact that when the volume fraction occupied by 
the ionized gas becomes close to unity most photons are 
preferentially used to sustain the reached ionization level 
rather than to create new ionized regions. 

In principle, a higher photon injection in the IGM could 
result both in an increase (as larger HII regions are pro- 
duced) or a decrease (as the number of sources is reduced 
by the effect of radiative feedback) of the filling factor. 
Along the sequence run B, D, A, C the number of ioniz- 
ing photons injected in the IGM is progressively increased. 
Then Fig. 9a allows us to conclude that the former effect 
is dominating, i.e. the radiative feedback is of minor im- 
portance. 

For the specific case of run A (we do not discuss in 

further detail He reionization in this paper) we show the 
analogous filling factor for the doubly ionized He; the low- 
est curve in Fig. 9a shows its evolution. As Rs,He — O.li?^ 
(see eqs. [5] and [7]), the helium filling factor should be 
approximately 0.1 % of the hydrogen filling factor, but it 
increases when the effects of the ionized hydrogen sphere 
overlapping is included (see Fig. 9a) . The stellar spectrum 
contributed by the reionizing galaxies is relatively soft and 
the number of photons above the He"*""*" ionization thresh- 
old very limited. Hence, He reionization does not occur 
in our model and could be possibly achieved only through 
quasar-like sources with harder spectra. 

In Fig. 9b we analyze the effect of negative feedback and 
the normalization of the dark matter fluctuation power 
spectrum, erg. By comparing run A with run Al, in which 
only objects with masses M > Mh are considered, we 
see that, apart from very high redshifts where few such 
objects are present and the bulk of photon production is 
due to objects with M < Mh, the main contribution to 
the IGM reionization comes from sources with M > Mh 
and Pop His (with M < Mh) alone would not be able 
the reionize the IGM. This statement remains valid even 
in the absence of radiative feedback effects on Pop III as 
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concluded from the results of a run A in which radiative 
feedback has not been allowed. In this case, although a 
slightly higher filling factor is obtained as the contribu- 
tion from previously suppressed Pop III is now present, 
the difference with run A is never higher than 10%. 

Finally, in run A2, the normalization of the fluctuation 
spectrum has been decreased to = 0.5. The evolution 
of the filling factor appears similar to that in run A, but 
shifted towards lower redshift. This is due to the delayed 
evolution of this model with respect to the reference one, 
which results in a lower number of luminous objects at a 
given redshift. 

7.2. SUVB 

The derived evolution of the SUVB is presented in 
Fig. 10a at the mean frequency of the LW band, /ii/o=12.45 
eV, for runs C, A, D and B (from the top to the bottom). 
The intensity scales with the product ft x x fesc as de- 
duced from eqs. (2), (15) and (16). Typically, a flux lower 
than lO"^** erg cm^^ s~^ Hz~^ sr~^ is found at z ^ 15, 
except for case C for which it is found to be about one 
order of magnitude larger. As an intensity of at least few 
xlO^^^ erg cm^^ s~^ Hz"'^ sr~^ at the Lyman limit is 
needed for an appreciable negative feedback at these red- 
shifts (see Fig. 5), we conclude, in agreement with CPA, 
that the SUVB flux is too weak to affect small mass struc- 
ture formation at high redshift. At these early epochs, neg- 
ative feedback is instead induced by the direct flux from 
pre-existing nearby objects. The SUVB becomes intense 
and dominates the direct flux for 2;<15 (earlier for run C) 
and consequently it governs the formation of late forming 
Pop Ills. The additional cases for runs Al and A2 are 
shown in Fig. 10b. Run Al produces the same SUVB as 
run A at low redshift; at early times, where there is little 
contribution from objects with M > Mh, the SUVB in- 
tensity is lower. The case with a lower normalization (run 
A2) also produces a lower SUVB, as expected. 

7.3. SFR 

We have calculated the star formation rate per comoving 
volume (SFR) predicted by our model and we have com- 
pared it with the values derived by the most recent studies 
at lower redshift (Lilly et al. 1996; Connolly et al. 1997; 
Madau, Pozzetti & Dickinson 1998; Steidel et al. 1998). 
Inside each simulation box we have derived the SFR as 
SFR = J2i ^/''i) where V is the comoving volume of 
the simulation box and s/fj are the star formation rates 
of the individual luminous objects: 



sfn = 



(28) 



in units of Mq yr"-'^. Here Mi,^i is given by eq. (1), show- 
ing that the SFR crucially depends on the parameters 
and /*; is defined as the free-fall time for objects with 
mass M > Mhy, while for small mass objects (M < M^y), 
witnessing only a single burst of star formation, we assume 
it equal to the Hubble time: 



tn M < Mt 
tff M > Mb- 



(29) 



where tu = {2/3)Hq ^ (1 + z)~^/^. As we have already 
pointed out, the simulations, although resolving small 



mass objects, might miss halos larger than Mmax.s] this 
value increases with decreasing redshift. To take into ac- 
count this limitation and correct for it, we include also 
objects with masses M > Mmax,s and model their con- 
tribution as follows. At each redshift, this contribution is 
given by: 



sfrps 



te{M) 



Af{M)dM, 



(30) 



where Mi, is given in eq. (1), while M{M) and M„iax are 
defined in eq. (16). We implicitly assume that these ob- 
jects do not suffer from radiative feedback (this assump- 
tion is justified by the fact that these are usually objects 
with masses M > Mh) and that these halos host lumi- 
nous objects. We then add this contribution to the SFR 
derived from the simulations. The results are shown in 
Fig. 11a for runs A-D. Obviously, the higher the fb x /* 
product, the higher SFR is obtained. Actually, although 
run A and D use the same values for the above parame- 
ters, in run D more luminous objects are formed (as the 
SUVB is lower [see § 7.2] more objects escape negative 
feedback) and this results in a slightly higher SFR. We 
compare the curves obtained with the most recent mea- 
surements of the SFR, namely: Lilly et al. (1996) [circles], 
Madau, Pozzetti & Dickinson (1998) [triangles], Connolly 
et al. (1997) [squares] and Steidel et al. (1998) [crosses]. 
The recent lower limit to the SFR at zk2, set by SCUBA 
(Hughes et al. 1998) coincides with the point correspond- 
ing to Steidel et al. (1998). The points from Madau, 
Pozzetti & Dickinson (1998) are derived from a study of 
the Lyman break galaxies in the Hubble Deep Field (HDF) 
and show a decline of the SFR at high redshift. The mea- 
surements by Steidel et al. (1998) do not show the same 
decline and favour a more or less constant SFR at high 
redshift. This determinations are based on a study of star- 
forming galaxies at z >3.8 in a field much wider than the 
HDF, although shallower. All points are corrected for dust 
extinction as in Steidel et al. (1998). From Fig. 11a, we see 
that runs B and C can be taken respectively as lower and 
upper limits to the SFR obtained from our model, thus 
constraining the product fb x /^.. Run A and D produce 
a trend which apparently well matches the observations, 
suggesting a likely value of the product fb x fi, around 
psO.Ol. In Fig lib we study the effect of negative feedback 
and (Tg- From a comparison between run A and Al, where 
only objects with M > Mh are considered, we see that, 
while at high redshift the main contribution to the SFR 
comes from small mass objects (M < Mh), at lower red- 
shift, their contribution becomes negligible. Finally, run 
A2, with CT8=0.5, produces a lower SFR as less luminous 
objects are present at a given redshift, again the same 
point made in § 7.1. 

Additionally, we have derived the evolution of the stel- 
lar density parameter, fli,{z) = Pi,,tot{z)/ pent, produced 
by run A; pi,^tot{z) is the total mass density in stars at 
redshift z inside the comoving volume of the box V=12b 
Mpc^, and is given by 



P*,tot{t) 



SFR{t)dt, 



(31) 
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where ton is the Hubble time corresponding to Zori=29.6. 
We find that by the time the IGM is completely reion- 
ized, ^~10, only 2% of the stars observed at z = 
(O*(0) w 4.9 X 10-3 [Fukugita, Hogan & Peebles 1998]) has 
been produced. This apparently small amount of stars is 
sufficient to produce the necessary ionizing flux. The num- 
ber of ionizing photons produced by zwlO can be derived 
as follows. At ^«10, M^^tot « 3.8 x 10^ Mq, corresponding 
to an ionizing photon rate Si{0) ~ 2.62 x 10^^ s"^ for the 
standard case A (see eq. [4]). The time at which the lowest 
mass OB stars expire is w40 Myr (Oey & Clarke 1997), 
thus, the total number of ionizing photons escaping from 
luminous objects is Nj « 3.3 x 10^®. As the number of 
baryons in the box is Nf, = V^hPcrit/ [p-'mp) « 4.7 x 10^^, 
where /i « 1.27 is the mean molecular weight, there are 
Ni/N(,«7 ionizing photons per baryon available. At 2;«10 
the Hubble time is Ir ~ 3.59 x 10^ yr, while the mean hy- 
drogen recombination time over the time interval 2:«30-10 
of interest here is {tree) ~ 1.33 x 10* yr. Thus, typically an 
atom has recombined about tu / {tree) ~3 times. However, 
the ionizing photon budget provided by the formed stars is 
sufficient to balance this recombination and keep the gas 
ionized. 

A rough estimate of the mean metallicity produced by 
these stars can be calculated. Given the adopted IMF 
and the total mass in stars at 2;«10, M^^^ot, we can de- 
rive the number of SNe in the simulation box, and, as- 
suming that each SN produces wl Mq of heavy elements, 
we find that the total mass of heavy elements at zwlO is 
Mmet ~ 2.81 X 10^ Mq. As the total mass in baryons in 
the considered cosmic volume is Mb^tot ~ 5.06 x 10^^ Mq, 
we find that the mean metallicity at zwlO is: 



{Z) = 



Mr, 



Mb, 



5.5 X IQ-'^ « 3 X IQ-^Zq. 



(32) 



7.4. Galaxy evolutionary tracks 



In this Section we discuss the final fates of primordial 
objects and we show their relative numbers in Fig. 12. 
The straight lines represent, from the top to the bottom, 
the number of dark matter halos, dark objects, naked stel- 
lar clusters and gaseous galaxies, respectively. The dotted 
curve represents the number of luminous objects with large 
enough mass (M > Mh) to make the H line cooling effi- 
cient and become insensitive to the negative feedback. We 
remind that the naked stellar clusters are the luminous 
objects with M < M^y, while the gaseous galaxies are the 
ones with M > Mby-, thus, the number of luminous objects 
present at a certain redshift is given by the sum of naked 
stellar clusters and gaseous galaxies. We first notice that 
the majority of the huninoiis objects that are able to form 
at high redshift will experience blowaway, becoming naked 
stellar clusters, while only a minor fraction, and only at 
z ^ 15, when larger objects start to form, will survive and 
become gaseous galaxies. An always increasing mimber of 
luminous objects is forming with decreasing redshift. until 
2;~15, where a fiattening is seen. This is due to the fact 
that the dark matter halo mass function is still dominated 
by small mass objects, but a large fraction of them can- 
not form due to the following combined effects: i) towards 
lower redshift the critical mass for the collapse (Merit) in- 
creases and fewer objects satisfy the condition M > Merit', 
ii) the radiative feedback due to either the direct dissoci- 



ating flux or the SUVB (see § 4) increases at low redshift 
as the SUVB intensity reaches values significant for the 
negative feedback effect. When the number of luminous 
objects becomes dominated by objects with M > Mh, by 
zwlO the population of luminous objects grows again, ba- 
sically because their formation is now unaffected by neg- 
ative feedback. A steadily increasing number of objects 
is prevented from forming stars and remains dark; this 
population is about w99% of the total population of dark 
matter halos at zw8. This is also due to the combined ef- 
fect of points i) and ii) mentioned above. This population 
of halos which have failed to produce stars could be identi- 
fied with the low mass tail distribution of the dark galaxies 
that reveal their presence through gravitational lensing of 
quasars (Hawkins 1997; Jimenez et al. 1997). It has been 
argued that this population of dark galaxies outnumbers 
normal galaxies by a substantial amount, and Fig. 12 sup- 
ports this view. At the same time, CDM models predict 
that a large number of satellites (a factor of 100 more than 
observed [Klypin et al. 1999; Moore et al. 1999b]) should 
be present around normal galaxies. Many of them would 
form at redshifts higher than five and would survive merg- 
ing and tidal stripping inside larger halos to the present 
time. Their existence is then linked to that of small mass 
primordial objects and the natural question arises if these 
objects can be reconciled with the internal properties of 
halos of present day-galaxies. 

A question that naturally arises is which is the fate of 
the naked stellar objects. The actual number of naked 
stellar objects should not be much higher than the one 
at z f«8, for two reasons: (i) with decreasing redshift, a 
lower number of luminous objects is subject to blowaway, 
as increasingly larger objects form; (ii) as already pointed 
out, at redshifts below «6, our estimate of Mby becomes 
less accurate, actually larger than the real value (see FT), 
thus, the number of luminous objects undergoing a blow- 
away is decreased also by this effect. With this caveat, 
approximately 2 naked stellar objects should end up in 
the Galactic halo, where, at present time, only stars with 
masses in the range 0.1-1 M© are still shining. As the 
maximum total mass of a naked stellar object at z «8 is 
« 10^ Mq, given a Salpeter IMF, we find that the upper 
limit for the number of the above relic stars is « 4 x 10^. 
As in the Galaxy, given the same IMF, w 2 x 10^^ stars 
with masses in the range 0.1-1 Mq are present, one out 
of w 5 x 10** stars comes from naked stellar objects. The 
upper limit on the metallicity of such stars is the IGM one 
at z wlO, i.e. Z fa 3 x 10~'^Zq, as derived in eq. (32). 

We have compared the results for runs A and C. The 
only significant difference is that in case C, where a higher 
SUVB is found, the redshift at which the M > Mh pop- 
ulation peaks is higher, as from z « 15 there are very 
few objects that can avoid negative feedback and the only 
surviving halos are the ones for which H line cooling is 
efficient. 

8. COMPARISON WITH PREVIOUS WORK 

In the last decade there have been a number of different 

works dedicated to reionization. The "first generation" of 
studies used very simplified descriptions of the problem in 
terms of the treatment of the physical processes, cosmo- 
logical evolution of the ionizing sources, or both. These 
pioneering works (Fukugita & Kawasaki 1994; Tegmark, 
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Silk & Blanchard 1994; Giroux & Shapiro 1996) found 
results often in disagreement among them, probably due 
to the different effects/cosmological models included or ne- 
glected. More recent works have tackled again the problem 
improving on at least two crucial ingredients: a detailed 
cosmological hydrodynamic evolution of the IGM/galaxy 
formation, which gives a more detailed distribution of ion- 
izing sources and of the gas, and the inclusion of newly 
discovered feedback effects, along the lines presented here. 
GO performed detailed simulations for a CDM+A cosmo- 
logical model with an approximate treatment of the radia- 
tive transfer and allowing for IGM clumping. The N-body 
simulations used here have a larger number of particles 
(256^ instead of their 64/128^) although the size of the box 
is practically the same. Their spatial resolution (wl kpc) 
can only marginally resolve the HII spheres which, for Pop 
III objects, are often smaller than that size. Also, we treat 
in more detail the effects of radiative and stellar feedback, 
even if we cannot compete with the wealth of information 
on the IGM density distribution that can be extracted 
from their simulations. They conclude that the reioniza- 
tion redshift is ~7, i.e. about 3 redshift units below the 
one found here. The discrepancy might depend both on 
the different cosmological model adopted and/or on the 
above differences. A low value for the reionization red- 
shift is found also in the semi-analytical model presented 
by VS: for both an open and a critical density universe, 
they use a mass hmction alternative to the usually adopted 
PS formalism and, given a prescription for galaxy forma- 
tion and photon production, they find that the reionization 
is complete by redshift 6.8 for the open universe and 5.6 
for the critical density one. As in their calculation they 
include the treatment of the gas dumpiness, this could 
again explain the discrepancy with the value found in the 
present work. Additionally, deriving the radius of the HII 
region surrounding a source, they neglect the influence of 
the expansion of the universe, underestimating the real 
radius and obtaining a late reionization epoch. Haiman 
& Loeb (1997) used a semi-analytical model based on PS 
formalism which adopted a strong version of the radiative 
feedback (but not considering the stellar one), assuming 
the suppression of star formation inside objects with virial 
temperature below 10^ K. They also neglect the effect of 
gas dumpiness. This approach has the advantage of being 
very flexible when exploring the dependence of reionization 
on various parameters, but clearly lacks the crucial infor- 
mation about the spatial distribution of ionizing sources. 
They find that for a large range of CDM models reioniza- 
tion occurs at a redshift ^10, consistent with the value 
we find in the present work. The basic agreement with 
their result derives from the fact that we also find that 
reionization is mostly driven by objects collapsed through 
H line cooling. Stellar feedback instead introduces a differ- 
ent prediction for the fate of ionizing objects as discussed 
above. 

9. DISCUSSION 

In this paper we have studied the reionization of the 
universe due to an inhomogeneous distribution of sources, 
including a number of relevant physical processes and feed- 
back effects whose importance for this type of studies has 
only recently been recognized. Our approach provides a 
reliable picture of the actual process of IGM reionization. 



Maybe, the major feature not considered here is con- 
stituted by the possible IGM density inhomogeneities, 
as throughout the paper we have assumed an homoge- 
neous gas distribution. The IGM dumpiness can have 
in principle different effects on the reionization of the 
universe. Both GO and VS have studied the process 
of reionization taking into account the effect of the in- 
homogeneous gas distribution through a clumping factor 
^ — {p]gm) KpigmY 1 that, included in the calculation of 
the recombination process, increases the number of ioniz- 
ing photons necessary to reionize the IGM by the amount 
C, thus delaying the IGM complete reionization with re- 
spect to the homogeneous case. The gas dumpiness in- 
fluences also the shape of the ionized regions: indeed, 
the propagation of the ionizing front in an inhomogeneous 
medium, does not result in a spherically symmetric HII re- 
gion, with ionization proceeding in an inside-out fashion, 
but rather the reionization occurs outside-in, starting in 
voids and gradually penetrating in high-density regions as 
pointed out by Miralda-Escude. Haehnelt & Rees (1998). 
If this is the case, recombinations from overdense regions 
which are presumably very compact, would contribute to 
the average recombination rate only at later times, when 
most of the volume has already been ionized. Note that 
this effect is opposite to the first one, indicating that at 
a given redshift in the homogeneous case one might be 
underestimating the filling factor of the ionized gas. An 
even more refined approach to the problem should include 
the combined treatment of the complete radiative trans- 
fer equation in an inhomogeneous medium. The very first 
steps in this directions have been presented by Norman, 
Paschos & Abel (1998), who are trying to incorporate ra- 
diative transfer into 3D hydrodynamic cosmological sim- 
ulations. As a preliminary study of the effect of the gas 
dumpiness on our results, we have made additional runs 
in which we include a clumping factor in the calculation 
of the recombination process. For C we have adopted the 
curves given in GO and VS. In Fig. 13 the HII filling fac- 
tor for run A is shown, together with the curves derived 
in the presence of clumping. The GO clumping factor 
(dashed curve) produces little difference, as it is close to 
unity until z «10. On the other hand, the VS one (dotted 
line), reaching much higher values, results in a substan- 
tially lower filling factor and a delay in the IGM reioniza- 
tion, which is not complete by 2w8. 

Another important output of the model is the star for- 
mation history. A direct comparison with most of the 
previous calculations of the SFR, both in numerical and 
semi-analytical approaches, can only be approximate, as 
those mainly concentrated on lower redshift range than the 
ones studied here (see for example Guiderdoni et al. 1998; 
Baugh et al. 1998b; Pei, Fall & Hauser 1998; Nagamine, 
Gen & Ostriker 1999). VS, on the other hand, calculate 
the SFR up to 2;«20, obtaining values slightly lower both 
than ours and the ones observationally derived. 

Finally, the result that a high fraction of dark objects 
is present, is intriguing, as the question of whether dark 
galaxies can exist is a longstanding one (Dekel & Silk 
1986). Indeed, the fact that the luminosity function has a 
flatter faint slope than the associated halo mass function 
can be explained with an increasing fraction of dark galax- 
ies towards lower masses. Some authors (Babul & Rees 
1992; Metcalfe et al. 1997) claim that the initial burst of 
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the star formation in small mass objects may provide a 
population of disappearing dwarfs proposed for interpret- 
ing the faint galaxy counts. On the other hand, Jimenez 
et al. (1997) propose that a dark galaxy forms as a low- 
density disk in a dark halo of high spin parameter: such 
a galaxy may have a surface density too low to produce a 
significant star formation rate. Our model offers an alter- 
native explanation to the nature of such objects. 

A comparison of the predictions of our model is cur- 
rently possible only in a indirect way, as done through- 
out the paper for several quantities. However, in a near 
future, thanks to the availability of powerful spatial and 
ground-based instrumentation it will become possible to 
directly probe most of the activity occurring prior to com- 
plete reionization. For example, NGST is expected to pro- 
duce the detection of a large number of high-z SNe (Marri 
& Ferrara 1998) that could be used to trace early star for- 
mation and stellar feedback in Pop III objects, and possi- 
bly the IGM, via absorption experiments. NGST will be 
also able to determine with great accuracy the reionization 
epoch from the spectra of high-redshift sources (Haiman & 
Loeb 1998c). The ionized regions can be mapped directly 
in free-free emission (Oh 1999) by the Square Kilometer 
Array (SKA). Reionization should leave a measurable im- 
print on the CMB, as secondary anisotropics are produced 
via the kinetic Sunyaev-Zeldovich effect by inhomogeneous 
reionization (Knox, Scoccimarro & Dodelson 1998; Gruzi- 
nov & Hu 1998); this experiment is well into the capabili- 
ties of MAP and Planck missions, submillimetre telescopes 
will be used to search for dust, produced by primordial ob- 
jects, and possibly molecules; the dust can be an important 
reservoir of the metals synthesized by the first SNe and it 
might have important consequences for the exploration of 
the early universe. Finally, radio surveys might reveal sig- 
natures of reionization and give insight into the thermal 
evolution of neutral hydrogen through its redshifted HI 21- 
cm line emission. These features are within reach of new 
generation radio telescope like SKA (Madau, Meiksin & 
Rees 1997; Tozzi et al. 1999). Clearly, exciting times are 
ahead of us. 

10. SUMMARY 

We have studied inhomogeneous reionization in a crit- 
ical density CDM model, including a detailed treatment 
of radiative and stellar feedback processes on galaxy for- 
mation. The main results discussed in this paper can be 



summarized as follows. 

• Galaxies are able to reionize the neutral atomic hy- 
drogen by a redshift zwlO, while molecular hydro- 
gen is completely dissociated at very high redshift 
(zw25). 

• IGM reionization is basically driven by objects col- 
lapsed through H line cooling (M > Mh), while 
small mass objects (M < Mh) play only a minor 
role and even in the absence of a radiative nega- 
tive feedback they would not be able the reionize 
the IGM. 

• The soft-UV background intensity is too low to pro- 
duce sensible negative feedback effects on low mass 
galaxy formation at redshift z ^15, where the ra- 
diative feedback is dominated by the direct flux from 
pre-existing objects. At lower redshifts the SUVB is 
instead dominant and reaches interesting values to 
influence the subsequent galaxy formation process. 

• The evolution of the star formation rate obtained 

shows a trend consistent with the most recent mea- 
surements and the match constrains the baryon con- 
version efficiency into stars in a narrow range around 
psO.Ol. Only about 2% of the stars observed at = 
is required to reionize the universe. This corresponds 
to an average IGM metallicity at redshift zp^IO equal 
to (Z) w 3 X IQ-'^Zq. 

• A consistent fraction of halos is prevented from form- 
ing stars by either the condition M < Merit or the 
effect of radiative feedback; this population of dark 
objects reaches w99 % of the dark matter halo pop- 
ulation at ^;~8. 
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RUN 


h 


/* 


fesc 


cooling 


<^8 


A 


0.08 


0.15 


0.2 


H+Ha 


0.6 


B 


0.08 


0.05 


0.2 


H+Ha 


0.6 


C 


1.00 


0.15 


0.2 


H+Ha 


0.6 


D 


0.08 


0.15 


0.1 


H+Ha 


0.6 


Al 


0.08 


0.15 


0.2 


H 


0.6 


A2 


0.08 


0.15 


0.2 


H+Ha 


0.5 



Table 1 

Parameters of the calculation: fraction of virialized baryons that are able to cool and be available for star 
formation, /;,; star formation efficiency, /*; photon escape fraction, fesc; object cooling mechanism; dark 

matter fluctuation power spectrum normalization, as. 
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Fig. 1. — Evolution of the adopted emission spectrum as a function of photon energy, at four different evolutionary times t = (0, 1, 2, 32) X 
10^ yr. The specific luminosity corresponds to one solar mass of stars formed. 
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Fig. 2. — Evolution of the ionized atomic hydrogen fraction as a function of depth, for a Pop III at 2: = 29.6 and an incident flux 
Js,o = 1^21 IQ-^^ erg s^^ cm~-^ Hz~^ sr~^. The curves are derived for the BC steUar spectrum calculated at t = 2 X 10^ yr and, for 
comparison, for a power law spectrum (PL) with a=l.5. Different values of the /3 parameter are shown. The horizontal line represents the 
initial condition. 
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Fig. 3. — As Fig. 2 for the gas temperature. 
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Fig. 4. — As Fig. 2 for the neutral molecular hydrogen fraction. 
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Fig. 5. — Minimum total mass for self-shielding from an external incident flux with intensity J^^o a^t the Lyman limit. The curves are for 
different redshift: from the top to the bottom z =8.3, 11.4, 15.4, 18.0, 22.1, 29.6. Circles show the value of (open) and M^rit (filled). 
Radiative feedback works only in the solid portions of the curves. 
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Fig. 6. — Minimum total mass for collapse of an object as a function of redshift according to Tegmark et al. (1997) {Merit, solid line) and 
maximum mass for blowaway as a function of redshift {Mi,y, dotted line). Mi,y is shown for the three runs A, B, C (see Tab. 1). 
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Fig. 7. — Possible evolutionary tracks of objects as determined by the processes and feedbacks included in the model. Explanations are 
given in §6. 
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The figures are available on 
http:/ /www. arcetri.astro.it/ ~ferrara 



Fig. 8. — Global spatial structure of photoionized regions at 2=22.1 (a), 19.8 (b), 18.0 (c) and 15.4 (d). 
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Fig. 9. — (a) dissociated molecular hydrogen (upper set of lines), ionized atomic hydrogen (middle set) and doubly ionized helium (bottom 
line) filling factor as a function of redshift for different runs: A (solid line), B (dotted), C (dashed-dotted) and D (dashed), (b) same as (a) 
for runs: A (solid lino), Al (dotted) and A2 (dashed-dotted). 
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Fig. 10. — (a) SUVB intensity as a function of redshift, at the central frequency of the Lyman- Werner band hvQ=12.45 eV for runs A-D. 
(b) same as (a) for runs A, Al, A2. 
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Fig. 11. — (a) comoving star formation rate as a function of redshift, for runs C, D, A, B, from the top to the bottom, respectively. Data 
points are from Lilly et al. (1996) [circles]; Connolly et al. (1997) [squares]; Madau, Pozzetti Sz Dickinson (1998) [triangles] and Steidel et 
al. (1998) [crosses], (b) same as (a) for runs A, Al, A2. 
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Fig. 12. — (a) Number evolution of different objects in the simulation box for run A (see discussion in § 7.4) (b) same as (a) for run C. 
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Fig. 13. — Ionized atomic hydrogen filling factor as a function of redshift for diff'ercnt runs: A (solid line), A with the clumping factor curve 
of GO (dashed) and VS (dotted). 



